#ifndef __MESH_H__
#define __MESH_H__
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include "parmetis.h"
#include "constants.h"
#include "OpenPFEM.h"

#define PARDITION_NUM 20

typedef struct PTW_
{
    INT NumRows;     // 表示需要表示的元素的个数
    INT NumCols;     // 表示与元素相联系对象的总个数
    INT *Ptw;        // 长度为NumRows+1
    INT *Entries;    // 长度由非零元素来决定，可以理解成矩阵中的列号
    INT *AndEntries; // 第二个元素数组, 其位置定义与Entries相同
} PTW;

// 下面开始来进行具体的网格元素的定义，我们采用先定义几何元素，然后再来组合成网格元素的方式
/* 点的定义 */
typedef struct VERTEX_
{
    INT BD_ID;       // 边界信息
    INT Index;       // 在当前进程的局部编号
    DOUBLE *Coord; // 物理坐标
} VERT;
/* 线的定义 */
typedef struct LINE_
{
    INT Vert4Line[2]; // 两个顶点
    INT BD_ID;        // 边界信息
    INT Index;        // 在当前进程的局部编号
} LINE;
// 需要的函数：计算长度，切线，计算线上指定比例某点的坐标
/* 面的定义 */
typedef struct FACE_
{
    INT NumVerts;                // 点的个数
    INT NumLines;                // 边的个数
    INT Vert4Face[MAXVERT4FACE]; // 顶点
    INT Line4Face[MAXLINE4FACE]; // 边
    INT BD_ID;                   // 边界信息
    INT Index;                   // 在当前进程的局部编号
} FACE;
/* 体的定义 */
typedef struct VOLU_
{
    INT NumVerts;                // 点的个数
    INT NumLines;                // 边的个数
    INT NumFaces;                // 面的个数
    INT Vert4Volu[MAXVERT4VOLU]; // 顶点
    INT Line4Volu[MAXLINE4VOLU]; // 边
    INT Face4Volu[MAXFACE4VOLU]; // 面
    INT Index;                   // 在当前进程的局部编号
} VOLU;

typedef struct SHAREDGEO_
{
    INT SharedNum;    // 与本进程共享的节点个数
    INT *Index;       // 这些点在本进程中的局部编号
    INT *SharedIndex; // 这些点在共享进程中的编号
    INT *Owner;       // 每一个共享的几何元素的宿主进程号
} SHAREDGEO;

typedef struct SHAREDINFO_
{
    INT worlddim;
    MPI_Comm neighborcomm;
    INT NumNeighborRanks; // 与当前进程有共享关系的邻居进程号个数
    INT *NeighborRanks;   // 记录与当前进程有共享网格元素的邻居进程号
    // SharedVerts[i]: 与第i个邻居进程的共享节点信息
    SHAREDGEO *SharedVerts; // 存储网格中共享的节点信息
    SHAREDGEO *SharedLines; // 存储网格中共享的边的信息
    SHAREDGEO *SharedFaces; // 存储网格中共享的面的信息
    SHAREDGEO *SharedVolus; // 村粗网格中共享的体的信息(最多三维的情况下基本上用不到)
} SHAREDINFO;
/* 网格的定义 */
typedef struct MESH_
{
    INT worlddim;
    INT num_vert; // 点的个数
    INT num_line; // 边的个数
    INT num_face; // 面的个数
    INT num_volu; // 体的个数
    VERT *Verts;
    LINE *Lines;
    FACE *Faces;
    VOLU *Volus;
    INT *Fathers;   // 用来存储每个单元的父亲单元编号
    INT *Ancestors; // 用来存储每个单元的祖先单元

    //lhc 添加
    INT level_id;   // 用来存储本网格为第几层网格 用于判断网格之间的关系 祖先网格为第0层
    BOOL IsJacobiInfoAvailable; // 用来存储下面的三个信息是否都计算过 是为1 否为0
    DOUBLE *JacobiOfAllElem;    // 用来储存每个单元的Jacobi矩阵
    DOUBLE *InvJacobiOfAllElem; // 用来储存每个单元的Jacobi矩阵的逆矩阵
    DOUBLE *VolumnOfAllElem;    // 用来存储每个单元的体积

#if MPI_USE
    MPI_Comm comm;          // 本网格所在的通信域
    SHAREDINFO *SharedInfo; // 网格的分布式存储信息
#endif
} MESH;
//-------------------------------------------------------------------
//=====下面两个数据结构体目前只应用与网格划分过程中生成ShareInfo信息的时候使用
// 描述每个进程上某种共享元素的所有信息
typedef struct SHAREDINFO4GEO_
{
    INT SharedNum;    // 表示本进程中共享的几何元素个数
    INT *Owners;      // 表示本进程中共享几何元素所在的宿主进程，维数：SharedNum
    INT *Ptw;         // 表示每个共享元素所在的位置: 维数: SharedNum+1
    INT *SharedRanks; // 表示每个共享元素所在的进程: 维数: Ptw[SharedNum]
    INT *LocalIndex;  // 表示每个共享元素在所在进程的局部编号, 维数:Ptw[SharedNum]
} SHAREDINFO4GEO;
// 每个进程用来记录共享几何元素的结构体，其中中SharedRanks和LocalIndex是一一对应的

// 定义用来存储在每个进程上针对每个进程的共享信息结构体
typedef struct SHAREDGEO4RANK_
{
    INT SharedNumRanks; // 有多少个邻居进程与本进程共享几何元素
    INT *SharedRanks;   // 共享几何元素的进程号
    INT *Ptw;           // 针对每个进程在SharedIndex和Index中的位置
    INT *SharedIndex;   // 共享元素在相应进程中的局部编号
    INT *Index;         // 共享元素在本进程中的局部编号
    INT *Owners;        // 共享元素的宿主进程
} SHAREDGEO4RANK;
// 本结构体在每个进程中用来存储针对每个进程的共享信息, 便于信息的发送

// 6.25lhc修改， 与刘昊宸讨论一下 20221206
typedef struct
{
    DOUBLE level_sum[PARDITION_NUM + 1];
    INT level_num[PARDITION_NUM];
} LEVELSUMNUM;

// 针对一个DOUBLE一个INT组合的操作符 对应的操作均为求和
void dsumisumOP(LEVELSUMNUM *, LEVELSUMNUM *, int *, MPI_Datatype *);
// 针对DOUBLE[3]的运算操作符 每个位置分别对应运算 求和 求最大值 求最小值
void dsummaxminOP(DOUBLE *, DOUBLE *, int *, MPI_Datatype *);

/* creation */
void MeshCreate(MESH **mesh, INT worlddim, MPI_Comm comm);
void VertsCreate(VERT **verts, INT num, INT worlddim);
void LinesCreate(LINE **lines, INT num);
void FacesCreate(FACE **faces, INT num);
void VolusCreate(VOLU **volus, INT num);
void VertsReCreate(VERT **verts, INT num, INT new_num, INT worlddim);
void LinesReCreate(LINE **lines, INT num, INT new_num);
void FacesReCreate(FACE **faces, INT num, INT new_num);
void VolusReCreate(VOLU **volus, INT num, INT new_num);
void PTWCreate(PTW **ptw);
void SharedInfoCreate(MESH *mesh);
void SharedInfoUpdate(MESH *mesh, INT geotype, SHAREDINFO4GEO sharedgeo);
void SharedGeoCreate(SHAREDGEO **SharedGeo, INT num, INT *size);
void SharedGeoRenew(SHAREDGEO **SharedGeo, INT oldnum, INT newnum);

/* reading */
// 在root进程读入一个指定格式和网格类型的文件
// 下面的程序是用来读入Matlab格式的二维网格文件
// 需要介绍一下网格的格式！！！！
void Read_Data_line_root(MESH *mesh, const char *filename);
void Read_MatlabData_tri_root(MESH *mesh, const char *filename);
void Read_Data_tethedral_root(MESH *mesh, const char *filename);
// 网格信息的完备化
// 根据网格中的节点坐标、单元的节点编号、以及单元的三个表面的边界信息来生成一个完备的网格。
void Mesh3DComplete(MESH *mesh, INT *ElemBoundaries);
// 读入文件并建立一个网格对象
void MeshBuild(MESH *mesh, const char *filename, MESHDATATYPE meshdttp, MESHTYPE meshtype);
// 下面这个程序是将二维网格的每个三角形上的点和线进行重新排列，使得第一条边是主边
void Mesh2DLabel(MESH *mesh);
// 下面这个程序是将三维四面体网格上每个四面体上的点和线进行重排，使得第一条边和第四条边是最长边
void Mesh3DLabel(MESH *mesh);
// 将一个网格转变成可以进行Bisection的格式
void Mesh3D4BisectionLabel(MESH *mesh);
/* distribution */
void MeshPartition(MESH *mesh);
// 调用metis生成单元划分
void MetisPartition(MESH *mesh, PTW *Vert4Elems, INT NumElems, INT *ElemRanks);
// 得到没给进程上单元的几何
void GetElem4Ranks(INT *ElemRanks, INT NumElems, PTW *Elem4Ranks);
// GeoDim表示GEO元素的维数,0表示点，1表示线，2表示面，3表示体
void GetGEO4Elems(MESH *mesh, INT ElemDim, INT GeoDim, PTW *GEO4Elems);
/* output */
void Mesh2DPrint(MESH *mesh, INT printrank); // 在printrank进程打印网格信息
void Mesh3DPrint(MESH *mesh, INT printrank);
void MeshPrint(MESH *mesh); // 可以同时输出二维和三维的网格
void SharedMeshPrint(MESH *mesh, INT printrank);
void Mesh2DMatlabWrite(MESH *mesh, char *file);
// 按vtk格式输出网格信息
void Mesh3DVTKWrite(MESH *mesh, char *file);
/* refinement */
void MeshUniformRefine(MESH *mesh, INT time);
void MeshUniformRefine1D(MESH *mesh);
void MeshUniformRefine2D(MESH *mesh);
void MeshAdaptiveRefine(MESH *mesh, DOUBLE *ElemErrors, DOUBLE theta);
void MeshAdaptiveRefine1D(MESH *mesh, DOUBLE *ElemErrors, DOUBLE theta);
void MeshAdaptiveRefine2D(MESH *mesh, DOUBLE *ElemErrors, DOUBLE theta);
void MeshBisectionRefine2D(MESH *mesh, INT Num_Refine_Elems, INT *Refine_Elems);
void FindRefineElem(MESH *mesh, DOUBLE *PosterioriError, DOUBLE theta, INT *certain_refine_num, INT *REFINE_ELEMS);
void MeshUniformRefine3D(MESH *mesh);

// lhc 23.01.28 
// 将当前网格设置为祖先网格 使用该网格进行加密后的网格存储的ancestor为对应的该网格中的单元编号
void MeshSetAsAncestor(MESH *mesh);

/* free */
void MeshDestroy(MESH **mesh);
void SharedInfo4GEODestroy(SHAREDINFO4GEO **SharedInfo4GEO);
void SharedInfoDestroy(SHAREDINFO **SharedInfo);
void SharedGeoDestroy(SHAREDGEO *SharedGeo);

void PTWTranspose(PTW *GEO4Elems, PTW *GEO2Elems);
void GetGEO4Elems(MESH *mesh, INT ElemDim, INT GeoDim, PTW *GEO4Elems);
// 得到与每个几何元素相关联的单元集合
//  ElemDim表示单元的维数,1表示是一维网格，2表示是二维网格，3表示是三维网格
//  GeoDim表示GEO元素的维数,0表示点，1表示线，2表示面，3表示体
void GetGEO2Elems(MESH *mesh, INT ElemDim, INT GeoDim, PTW *GEO2Elems);
void GetGEO2Ranks(MESH *mesh, INT NumRanks, INT *ElemRanks, INT ElemDim, INT GeoDim, PTW *GEO2Ranks);
void GetGEO4Ranks(MESH *mesh, INT NumRanks, INT *ElemRanks, INT ElemDim, INT GeoDim, PTW GEO4Ranks);
void GetSharedINFO4GEO(PTW *GEO2Ranks, PTW *GEO4Ranks, SHAREDINFO4GEO **SharedINFO4GEO);
void ColletSharedGEO4Ranks(SHAREDINFO4GEO *SharedInfo4GEO, INT myrank, SHAREDGEO4RANK *SharedGEO4Ranks);
// 下面的函数是获得每个进程上每个单元上的几何元素在宿主进程中的编号
// GeoDim表示GEO元素的维数,0表示点，1表示线，2表示面，3表示体
void GetOwnerGEO4Elems(MESH *mesh, INT ElemDim, INT GeoDim, PTW *OwnerGEO4Elems);
// 下面的函数是获得每个单元在其宿主进程上所包含的局部几何元素的编号，即几何元素在单元进程上的编号
// GeoDim表示GEO元素的维数,0表示点，1表示线，2表示面，3表示体
void GetGEO4ElemsOnOwner(MESH *mesh, INT ElemDim, INT GeoDim, PTW *GEO4ElemsOnOwner);

void GetLocalGEO4Elems(PTW *GEO4Elems, PTW *Elem4Ranks, PTW *GEO2Ranks, PTW **LocalGEO4Elems);
void GetLine2Faces(MESH *mesh, PTW *PTWLine2Faces);
// 获得与每条线相联系的面的集合
void GetLine2Volus(MESH *mesh, PTW *PTWLine2Volus);
void PTWPrint(PTW *ptw);
void PTWPrint2(PTW *ptw);
void SharedInfo4GEOPrint(SHAREDINFO4GEO *SharedInfo4GEO);

INT MaxPosition(DOUBLE *VEC, INT start, INT end);
void Mesh3DCheck(MESH *mesh);
void BDIDCheck(MESH *mesh);
// 获得每两个点相对应的线的编号
void GetVert2Lines(MESH *mesh, PTW *Vert2Lines);
void PTWDestroy(PTW **ptw);

// 定义跟网格加密相关的数据结构
typedef struct REFINEMESH_
{
    // 存储网格
    MESH *Mesh;
    // 存储与每条边相联系的面的信息，稀疏矩阵的方式
    PTW *Line2Faces;
    // 存储每条边相联系的体的信息
    PTW *Line2Volus;
    // 记录每个单元上的每个局部面的主边在单元上的局部编号，长度4*num_volus
    INT *FaceLine4Volus; // 记录每个单元中的每个面的主边编号
    INT *ElemIndicator;  // 记录每个单元的加密信息
    INT *FaceIndicator;  // 记录每个面的加密信息：1表示加密，0表示不加密
    INT *LineIndicator;  // 记录每条边的加密信息
    INT *RefineElems;    // 记录需要加密的单元的编号
    INT *RefineFaces;    // 记录需要加密的面的编号
    INT *RefineLines;    // 记录需要加密的线的编号
    INT *FaceChildren;   // 记录每个面加密之后得到的子面集合(4*num_face)
    INT *FaceFather;     // 记录每个面的父亲面编号
    // 记录每个边的加密信息，知道当前边加密之后产生新边的编号
    // 第k条边的信息: [v[2k],v[2k+1]]分别表示第k条边加密之后得到的新的边的编号和相应的中点编号
    INT *RefineLineInformation;
    // 记录每个单元和面上每条边加密得到的中点编号，如果该边没有被加密就为-1
    INT *VoluLineMid;    // 数组的长度 6*num_volu
    INT *FaceLineMid;    // 数组的长度 3*num_face
    INT *FaceLineMidOld; // 记录四面体加密之前每个面上的三条边的中点编号
    // 存储网格中点到线的信息(可查知两个点所形成的线的编号)
    PTW *Vert2Lines;
    INT *Edge_Child_Vert; // 记录根据边产生的新点的编号
    INT *Edge_Child_Line; // 记录根据边产生的新非继承边的编号
    INT *Face_Child_Line; // 记录根据面产生的新线的编号
    INT *Face_Child_Face; // 记录根据面产生的新非继承面的编号
} REFINEMESH;
void MeshAdaptiveRefine3D(MESH *mesh, DOUBLE *ElemErrors, DOUBLE theta);
void GetFaceLine4Vlou(REFINEMESH *refinemesh);
void GetVoluLineMid(REFINEMESH *refinemesh, INT num_volu);
void BisectionRefineMesh(REFINEMESH *refinemesh, INT endvolu, INT endface, INT endline);
void GetRefineElements(DOUBLE *ElemErrors, DOUBLE theta, INT num_volu, INT *NumRefines, INT *RefineElems);
void GetFaceLineMid(REFINEMESH *refinemesh, INT Oldnum_face);
void BisectionLine(REFINEMESH *refinemesh, INT refineind, INT *newindices);
void BisectionTriangle(REFINEMESH *refinemesh, INT refineind, INT *newindices);
void BisectionTethedral(REFINEMESH *refinemesh, INT refineind, INT *newindices);
void CompleteAdaptiveMesh(REFINEMESH *refinemesh);
INT FindFacePosition(FACE face, INT Verts[], INT Positions[], INT num);
INT GetLocalPosition(FACE *Faces, INT *FaceChildren, INT father, INT Verts[], INT Positions[]);
void RefineMeshDestroy(REFINEMESH **refinemesh);
// 硬拷贝网格
MESH *MeshDuplicate(MESH *mesh);
SHAREDINFO *SharedInfoDuplicate(SHAREDINFO *sharedinfo);
//==================================================================
// 下面这些函数目前值被用于进行网格划分阶段的通讯功能，用户可忽视
//==================================================================
// 生成包裹的大小信息
INT PackageSizeVertInfo(MESH *mesh, PTW *Vert4Ranks, SHAREDINFO4GEO **SharedInfo4Verts, INT size,
                        INT **vertcount, INT **vertdisplc);
INT PackageSizeLineInfo(MESH *mesh, PTW *Line4Ranks, SHAREDINFO4GEO **SharedInfo4Lines, INT size,
                        INT **LineCount, INT **LineDisplc);
INT PackageSizeFaceInfo(MESH *mesh, PTW *Face4Ranks, SHAREDINFO4GEO **SharedInfo4Faces, INT size,
                        INT **FaceCount, INT **FaceDisplc);
INT PackageSizeVoluInfo(MESH *mesh, PTW *Volu4Ranks, INT size, INT **VoluCount, INT **VoluDisplc);
// 生成包裹
void PackageVertInfo(MESH *mesh, PTW *Vert4Ranks, SHAREDINFO4GEO **SharedInfo4Verts, INT size,
                     char **VertPackage, INT totalsize);
void PackageLineInfo(MESH *mesh, PTW *Line4Ranks, SHAREDINFO4GEO **SharedInfo4Lines, PTW **LocalVert4Lines,
                     INT size, char **LinePackageRoot, INT totalsize);
void PackageFaceInfo(MESH *mesh, PTW *Face4Ranks, SHAREDINFO4GEO **SharedInfo4Faces, PTW **LocalVert4Faces, PTW **LocalLine4Faces,
                     INT size, char **FacePackage, INT totalsize);
void PackageVoluInfo(MESH *mesh, PTW *Volu4Ranks, PTW **LocalVert4Volus, PTW **LocalLine4Volus, PTW **LocalFace4Volus,
                     INT size, char **VoluPackage, INT totalsize);
// 解包到网格
void UnpackedVert2Mesh(MESH *mesh, char *VertPackage, INT PackageSize);
void UnpackedLine2Mesh(MESH *mesh, char *LinePackage, INT PackageSize);
void UnpackedFace2Mesh(MESH *mesh, char *FacePackage, INT PackageSize);
void UnpackedVolu2Mesh(MESH *mesh, char *VoluPackage, INT PackageSize);
// 解包vert包裹中的shared信息
void UnpackedSharedInfo2Mesh(MESH *mesh, INT geotype, char *Package, INT PackageSize, INT *position);
// 产生线的全局编号
void LineGlobalIndexGenerate(MESH *mesh);
DOUBLE Distance4Verts(VERT vert0, VERT vert1, INT worlddim);
// 下面是根据加密单元来确定最终加密的线、面、体的集合
void FindRefineLineByElems(REFINEMESH *refinemesh, DOUBLE *ElemErrors, DOUBLE theta, INT *NumRefines);
void FindBisectionSetByRefineLines(REFINEMESH *refinemesh, INT *NumRefines);
// 寻找两个面的交线编号
void GetInterLine4Faces(FACE face0, FACE face1, INT *LineVerts);
// 寻找两个面的交点编号
INT GetInterVert4Lines(INT *Line0, INT *Line1);
// 寻找两个四面体的公共面的编号
INT GetInterFace4Volus(VOLU volu0, VOLU volu1);
// 根据边界类型调整网格低层级元素的边界ID
void MeshAdjustForDirichletBound(MESH *mesh, BOUNDARYTYPEFUNCTION *BoundType);
#endif